function [ber_uncertainty, capacity_loss, fig_handle] = channel_uncertainty(EbN0_dB, num_bits, uncertainty_levels)
% 信道不确定性影响分析
% 输入:
%   EbN0_dB - 比特信噪比范围 (dB)
%   num_bits - 仿真比特数
%   uncertainty_levels - 信道不确定性水平数组
% 输出:
%   ber_uncertainty - 误比特率矩阵
%   capacity_loss - 容量损失矩阵
%   fig_handle - 图形句柄

    EbN0_linear = 10.^(EbN0_dB/10);
    ber_uncertainty = zeros(length(EbN0_dB), length(uncertainty_levels));
    capacity_loss = zeros(length(EbN0_dB), length(uncertainty_levels));
    
    % 添加路径
    addpath('../Common');
    colors = color_definitions();
    
    fprintf('=== 信道不确定性影响分析 ===\n');
    
    for unc_idx = 1:length(uncertainty_levels)
        uncertainty = uncertainty_levels(unc_idx);
        fprintf('仿真信道不确定性 = %.2f...\n', uncertainty);
        
        for eb_idx = 1:length(EbN0_dB)
            fprintf('  Eb/N0 = %d dB...\n', EbN0_dB(eb_idx));
            
            % 生成随机比特
            tx_bits = randi([0, 1], num_bits, 1);
            
            % BPSK调制
            tx_symbols = 2 * tx_bits - 1;
            
            % 真实信道
            true_channel = sqrt(0.5) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
            
            % 估计信道 (含误差)
            channel_error = uncertainty * sqrt(0.5) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
            estimated_channel = true_channel + channel_error;
            
            % 添加AWGN噪声
            noise_power = 1 / (2 * EbN0_linear(eb_idx));
            noise = sqrt(noise_power) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
            
            % 接收信号
            rx_symbols = tx_symbols .* true_channel + noise;
            
            % 理想检测 (已知真实信道)
            ideal_symbols = rx_symbols ./ true_channel;
            ideal_bits = (real(ideal_symbols) > 0);
            ber_ideal = sum(tx_bits ~= ideal_bits) / num_bits;
            
            % 实际检测 (使用估计信道)
            detected_symbols = rx_symbols ./ estimated_channel;
            detected_bits = (real(detected_symbols) > 0);
            ber_uncertainty(eb_idx, unc_idx) = sum(tx_bits ~= detected_bits) / num_bits;
            
            % 计算容量损失
            % 理想容量
            ideal_snr = abs(true_channel).^2 * EbN0_linear(eb_idx);
            ideal_capacity = mean(log2(1 + ideal_snr));
            
            % 实际容量
            effective_snr = abs(estimated_channel).^2 * EbN0_linear(eb_idx);
            actual_capacity = mean(log2(1 + effective_snr));
            
            capacity_loss(eb_idx, unc_idx) = ideal_capacity - actual_capacity;
        end
    end
    
    % 绘制结果
    fig_handle = figure('Name', '信道不确定性影响', 'Position', [500, 500, 1200, 800]);
    
    % BER性能
    subplot(2,3,1);
    for unc_idx = 1:length(uncertainty_levels)
        semilogy(EbN0_dB, ber_uncertainty(:, unc_idx), ...
                 'Color', colors(unc_idx,:), 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('信道不确定性对BER的影响');
    legend(arrayfun(@(x) sprintf('σ=%.2f', x), uncertainty_levels, 'UniformOutput', false));
    
    % 容量损失
    subplot(2,3,2);
    for unc_idx = 1:length(uncertainty_levels)
        plot(EbN0_dB, capacity_loss(:, unc_idx), ...
                 'Color', colors(unc_idx,:), 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('容量损失 (bps/Hz)');
    title('信道不确定性导致的容量损失');
    legend(arrayfun(@(x) sprintf('σ=%.2f', x), uncertainty_levels, 'UniformOutput', false));
    
    % 信道估计误差分布
    subplot(2,3,3);
    uncertainty = 0.3;
    num_samples = 10000;
    true_chans = sqrt(0.5) * (randn(num_samples, 1) + 1i * randn(num_samples, 1));
    est_errors = uncertainty * sqrt(0.5) * (randn(num_samples, 1) + 1i * randn(num_samples, 1));
    
    % 计算误差角度和幅度
    error_angles = angle(est_errors ./ true_chans) * 180 / pi;
    error_magnitudes = abs(est_errors) ./ abs(true_chans);
    
    histogram(error_angles, 50, 'Normalization', 'pdf');
    hold on;
    
    % 理论分布 (假设高斯误差)
    x_theory = -180:1:180;
    angle_pdf = 1 / (uncertainty * sqrt(2*pi)) * exp(-x_theory.^2 / (2 * uncertainty^2 * 180^2 / pi^2));
    plot(x_theory, angle_pdf, 'r-', 'LineWidth', 2);
    
    grid on;
    xlabel('相位误差 (度)');
    ylabel('概率密度');
    title(sprintf('信道估计误差分布 (σ=%.2f)', uncertainty));
    legend('仿真直方图', '理论分布');
    
    % 不同检测算法性能
    subplot(2,3,4);
    detection_methods = {'理想检测', 'MMSE检测', 'ML检测', '鲁棒检测'};
    uncertainty = 0.2;
    ber_methods = zeros(length(EbN0_dB), length(detection_methods));
    
    for eb_idx = 1:length(EbN0_dB)
        tx_bits = randi([0, 1], num_bits, 1);
        tx_symbols = 2 * tx_bits - 1;
        
        % 真实和估计信道
        true_channel = sqrt(0.5) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
        channel_error = uncertainty * sqrt(0.5) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
        estimated_channel = true_channel + channel_error;
        
        noise_power = 1 / (2 * EbN0_linear(eb_idx));
        noise = sqrt(noise_power) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
        rx_symbols = tx_symbols .* true_channel + noise;
        
        % 1. 理想检测
        ideal_symbols = rx_symbols ./ true_channel;
        ideal_bits = (real(ideal_symbols) > 0);
        ber_methods(eb_idx, 1) = sum(tx_bits ~= ideal_bits) / num_bits;
        
        % 2. MMSE检测
        mmse_symbols = rx_symbols .* conj(estimated_channel) ./ (abs(estimated_channel).^2 + noise_power);
        mmse_bits = (real(mmse_symbols) > 0);
        ber_methods(eb_idx, 2) = sum(tx_bits ~= mmse_bits) / num_bits;
        
        % 3. ML检测 (简化)
        % 考虑信道不确定性，使用保守估计
        conservative_channel = estimated_channel * (1 - uncertainty);
        ml_symbols = rx_symbols ./ conservative_channel;
        ml_bits = (real(ml_symbols) > 0);
        ber_methods(eb_idx, 3) = sum(tx_bits ~= ml_bits) / num_bits;
        
        % 4. 鲁棒检测
        % 使用信道误差统计信息
        effective_noise = noise + tx_symbols .* (true_channel - estimated_channel);
        effective_noise_power = noise_power + uncertainty^2;
        robust_symbols = rx_symbols .* conj(estimated_channel) ./ (abs(estimated_channel).^2 + effective_noise_power);
        robust_bits = (real(robust_symbols) > 0);
        ber_methods(eb_idx, 4) = sum(tx_bits ~= robust_bits) / num_bits;
    end
    
    for method_idx = 1:length(detection_methods)
        semilogy(EbN0_dB, ber_methods(:, method_idx), ...
                 'Color', colors(method_idx,:), 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('不同检测算法性能 (σ=0.2)');
    legend(detection_methods);
    
    % 训练序列长度影响
    subplot(2,3,5);
    pilot_lengths = [8, 16, 32, 64];
    uncertainty = 0.3;
    ber_pilot = zeros(length(EbN0_dB), length(pilot_lengths));
    
    for pilot_idx = 1:length(pilot_lengths)
        pilot_length = pilot_lengths(pilot_idx);
        fprintf('仿真训练序列长度 = %d...\n', pilot_length);
        
        for eb_idx = 1:length(EbN0_dB)
            % 生成训练序列
            pilot_symbols = ones(pilot_length, 1);
            
            % 真实信道 (假设在训练期间保持不变)
            true_channel = sqrt(0.5) * (randn(1, 1) + 1i * randn(1, 1));
            
            % 接收训练序列
            noise_power = 1 / (2 * EbN0_linear(eb_idx));
            pilot_noise = sqrt(noise_power) * (randn(pilot_length, 1) + 1i * randn(pilot_length, 1));
            rx_pilots = pilot_symbols * true_channel + pilot_noise;
            
            % LS信道估计
            estimated_channel = mean(rx_pilots ./ pilot_symbols);
            
            % 计算估计误差
            estimation_error = abs(estimated_channel - true_channel);
            effective_uncertainty = uncertainty + estimation_error / abs(true_channel);
            
            % 数据检测
            data_bits = randi([0, 1], num_bits, 1);
            data_symbols = 2 * data_bits - 1;
            data_noise = sqrt(noise_power) * (randn(num_bits, 1) + 1i * randn(num_bits, 1));
            rx_data = data_symbols * true_channel + data_noise;
            
            detected_data = rx_data ./ estimated_channel;
            detected_bits = (real(detected_data) > 0);
            ber_pilot(eb_idx, pilot_idx) = sum(data_bits ~= detected_bits) / num_bits;
        end
    end
    
    for pilot_idx = 1:length(pilot_lengths)
        semilogy(EbN0_dB, ber_pilot(:, pilot_idx), ...
                 'Color', colors(pilot_idx,:), 'LineWidth', 2);
        hold on;
    end
    grid on;
    xlabel('Eb/N0 (dB)');
    ylabel('误比特率 (BER)');
    title('训练序列长度影响 (σ=0.3)');
    legend(arrayfun(@(x) sprintf('N_p=%d', x), pilot_lengths, 'UniformOutput', false));
    
    % 自适应调制
    subplot(2,3,6);
    modulation_schemes = {'BPSK', 'QPSK', '16QAM'};
    uncertainty_levels = [0.1, 0.2, 0.3];
    
    for mod_idx = 1:length(modulation_schemes)
        modulation = modulation_schemes{mod_idx};
        
        for unc_idx = 1:length(uncertainty_levels)
            uncertainty = uncertainty_levels(unc_idx);
            
            % 根据信道不确定性选择调制方式
            if uncertainty < 0.15
                bits_per_symbol = 1; % BPSK
            elseif uncertainty < 0.25
                bits_per_symbol = 2; % QPSK
            else
                bits_per_symbol = 4; % 16QAM
            end
            
            % 生成数据
            total_bits = num_bits * bits_per_symbol;
            tx_bits = randi([0, 1], total_bits, 1);
            
            % 调制
            if bits_per_symbol == 1 % BPSK
                tx_symbols = 2 * tx_bits - 1;
            elseif bits_per_symbol == 2 % QPSK
                tx_symbols = qpsk_modulate(tx_bits);
            else % 16QAM
                tx_symbols = qam16_modulate(tx_bits);
            end
            
            % 信道
            true_channel = sqrt(0.5) * (randn(length(tx_symbols), 1) + 1i * randn(length(tx_symbols), 1));
            channel_error = uncertainty * sqrt(0.5) * (randn(length(tx_symbols), 1) + 1i * randn(length(tx_symbols), 1));
            estimated_channel = true_channel + channel_error;
            
            noise_power = 1 / (2 * EbN0_linear(5) * bits_per_symbol); % 固定SNR点
            noise = sqrt(noise_power) * (randn(length(tx_symbols), 1) + 1i * randn(length(tx_symbols), 1));
            rx_symbols = tx_symbols .* true_channel + noise;
            
            % 检测
            detected_symbols = rx_symbols ./ estimated_channel;
            
            if bits_per_symbol == 1 % BPSK
                detected_bits = (real(detected_symbols) > 0);
            elseif bits_per_symbol == 2 % QPSK
                detected_bits = qpsk_demodulate(detected_symbols);
            else % 16QAM
                detected_bits = qam16_demodulate(detected_symbols);
            end
            
            ber_adaptive(unc_idx, mod_idx) = sum(tx_bits ~= detected_bits) / total_bits;
        end
    end
    
    % 绘制结果
    bar_data = reshape(ber_adaptive, [], 1);
    bar(bar_data);
    grid on;
    xlabel('调制方案');
    ylabel('误比特率 (BER)');
    title('自适应调制性能 (SNR=15 dB)');
    
    % 设置x轴标签
    x_labels = {};
    for mod_idx = 1:length(modulation_schemes)
        for unc_idx = 1:length(uncertainty_levels)
            x_labels{end+1} = sprintf('%s,σ=%.1f', modulation_schemes{mod_idx}, uncertainty_levels(unc_idx));
        end
    end
    set(gca, 'XTick', 1:length(bar_data), 'XTickLabel', x_labels);
    xtickangle(45);
    
    % 保存图片
    saveas(fig_handle, 'results/channel_uncertainty_analysis.png');
end

function symbols = qpsk_modulate(bits)
    symbols = zeros(length(bits)/2, 1);
    for i = 1:length(symbols)
        bit_pair = bits(2*i-1:2*i);
        if bit_pair == [0, 0]
            symbols(i) = 1 + 1i;
        elseif bit_pair == [0, 1]
            symbols(i) = -1 + 1i;
        elseif bit_pair == [1, 0]
            symbols(i) = 1 - 1i;
        else
            symbols(i) = -1 - 1i;
        end
    end
    symbols = symbols / sqrt(2);
end

function bits = qpsk_demodulate(symbols)
    bits = zeros(2*length(symbols), 1);
    for i = 1:length(symbols)
        s = symbols(i) * sqrt(2);
        if real(s) > 0 && imag(s) > 0
            bits(2*i-1:2*i) = [0; 0];
        elseif real(s) < 0 && imag(s) > 0
            bits(2*i-1:2*i) = [0; 1];
        elseif real(s) > 0 && imag(s) < 0
            bits(2*i-1:2*i) = [1; 0];
        else
            bits(2*i-1:2*i) = [1; 1];
        end
    end
end

function symbols = qam16_modulate(bits)
    num_symbols = length(bits)/4;
    symbols = zeros(num_symbols, 1);
    constellation = [-3-3i, -3-1i, -3+3i, -3+1i, ...
                       -1-3i, -1-1i, -1+3i, -1+1i, ...
                        3-3i, 3-1i, 3+3i, 3+1i, ...
                        1-3i, 1-1i, 1+3i, 1+1i] / sqrt(10);
    
    for i = 1:num_symbols
        % 获取4个比特
        bit_chunk = bits(4*i-3:4*i);
        % 转换为十进制索引 (0-15)
        symbol_idx = 0;
        for j = 1:4
            symbol_idx = symbol_idx + bit_chunk(j) * 2^(4-j);
        end
        % 映射到星座点 (1-16)
        symbols(i) = constellation(symbol_idx + 1);
    end
end

function bits = qam16_demodulate(symbols)
    constellation = [-3-3i, -3-1i, -3+3i, -3+1i, ...
                     -1-3i, -1-1i, -1+3i, -1+1i, ...
                      3-3i, 3-1i, 3+3i, 3+1i, ...
                      1-3i, 1-1i, 1+3i, 1+1i] / sqrt(10);
    
    bits = zeros(4*length(symbols), 1);
    for i = 1:length(symbols)
        % 找到最近的星座点
        distances = abs(symbols(i) - constellation);
        [~, closest_idx] = min(distances);
        
        % 转换为比特
        bit_pattern = de2bi(closest_idx - 1, 4)'; % 4比特
        bits(4*i-3:4*i) = bit_pattern;
    end
end